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ABSTRACT 



The analysis of stresses induced by contact between two bodies is inherently din'icult 
because the size of the contact zone is unknown and constantly changing throughout 
loading. To overcome these dilTiculties, two approximation methods have been devel- 
oped to determine the magnitude of contact stresses using the Rayleigh-Ritz method and 
the finite element method. Numerical optimization methods are employed to solve the 
contact problem. The solution techniques are compared to known analytical solutions 
and shown to yield accurate results. An application of this approach to solving the 
contact problem is illustrated by examining the response of a clamped sandwich com- 
posite beam to low velocity impact. It was found that the maximum shear stress is in- 
sensitive to lamina thickness, however an increase in the contact layer thickness resulted 
in a reduction in interfacial shear stress. In addition, it was noted that a nonlinear 
bending stress distribution in the contact layer intensified as the thickness of this layer 
increased. This phenomenon was found to be localized to the region of contact. Finally, 
it was found that the compressive transverse normal stresses increased as the thickness 
of the contact lamina increased. 
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I. INTRODUCTION 



A. MOTIVATION 

Contact stresses occur when two bodies exert forces over limited contact regions. 
Examples of contact stress include meshing gear teeth, cam shaft and pushrod contact, 
rollers in plate forming operations, roller and ball bearings in contact with races, shaft 
and journal bearing contact, and plate-pin connections. 

Contact zones can be point, line, or surface in geometry’. Because of the limited 
contact zone, the local stresses can be sufTiciently high to be of major concern to the 
designer. Consequently, a thorough understanding of this phenomenon is essential. The 
first successful analytical solution to the contact problem between two spheres was 
solved in the late 19th century by Hertz. His solution can only be applied to simple 
geometries such as spheres, cylinders, and flat plates. Because of these limitations, al- 
ternate solution techniques were needed to accommodate more complex geometries and 
boundary conditions. 

Unfortunately, the contact problem is difficult to study. The most significant diffi- 
culty is that the size of the contact zone is unknown and constantly changing throughout 
loading. Solutions are obtained by an iterative process. Consequently, the problem is 
highly nonlinear in its behavior. It is because of these difficulties that approximation 
and numerical methods are preferred in the solution of the contact problem. There are 
two general approaches used in the approximate techniques to solve the contact prob- 
lem. One class applies specific iterative procedures to solve a nonlinear system of 
equations that represent the contact state. For example, the contact condition can be 
simulated by the introduction of additional coupling terms into the system of equations. 
The other class constructs a functional that includes the contact body constraint. The 
functional is then minimized using specific numerical programming techniques. .Al- 
though these two classes use different procedures, there are several methods to incor- 
porate the contact condition into the problem formulation that are common to both 
classes. Examples of these include the penalty method and the augmented Lagrange 
multiplier method. These methods specify the manner in which the contact boundary 
conditions are treated. 

The objective of this study has two parts. First, two approximate solution tech- 
niques will be developed to obtain the solution of the general contact problem. These 



1 



techniques belong to the general class of methods that calculate a specific functional and 
then applies numerical programming methods to solve the contact problem. Second, 
these solution techniques will be verified by comparison with known analytical solutions. ) 
With the solution techniques verified, the methods will be applied to study an actual 
contact problem. 

I 

In order to accomplish the first objective, two models will be developed. The first | 
method utilizes the Rayleigh- Ritz method to solve the contact problem. .An assumed j 
deformation field that satisfies the boundarv' conditions is carefully selected. Theorv' of i 
elasticity relationships are applied to obtain an expression for the system's strain energy 
enabling calculation of the total potential. The minimization of the total potential en- 
ergy enables the calculation of the contact stresses at any point in the body. The second 
model developed uses finite element analysis to accomplish the same objectives. The 
numerical minimization technique used in both cases is the augmented Lagrange multi- 
plier method. 

To accomplish the second objective, both methods will be verified by comparison , 
with the Hertz solution of an infinitely long cylinder in contact with a flat plane. With 
verification completed, the contact stresses resulting from low velocity impact between 
objects and composite sandwich materials will be studied. 

B. LITERATURE SURVEY 

As discussed previously, there are two general approaches to solving the contact 
problem. One approach uses special procedures to solve a nonlinear system of 
equations. The other creates a functional that includes the contact boundarx' con- 
straints. Although different, both approaches use similar mathematical methods to in- 
corporate the contact condition into the problem formulation. Two of these methods 
are the Lagrange method and the penalty method. A detailed explanation of the method 
of employing the Lagrange method in nonlinear finite element analysis was discussed by 
Pian [Ref 1]. Alternately, the penalty finite element method was used by Cheng [Ref 
2] to solve the multibody contact problem. The fundamental concept of the latter ap- 
proach is the transformation of a constrained problem into an unconstrained one. This 
is done by penalizing (i.e., increasing) the objective function for constraint violations. 
.Although both methods are effective, there is substantial discussion on the limitations 
of both methods. Nour-Omid [Ref 3] described the positive and negative aspects of 
both methods. The Lagrange method has been shown to be the more accurate method. 
However, its usage requires the introduction of additional unknowns thus increasing the 



II. FORMULATION OF THE CONTACT PROBLEM 



A. PRINCIPLE OF MINIMUM TOTAL POTENTIAL ENERGY 

As discussed in the introduction, the limited utility of the analytical solutions ne- 
cessitated the development of solution techniques capable of handling the nonlinear be- 
havior of the contact problem with complicated geometr>' and complex boundary 
conditions. This study intends to develop two numerical procedures to solve the contact 
problem. In short, the procedures will use different methods to obtain a functional, the 
system's total potential energy, and then use sintilar methods to obtain the equilibrium 
condition. Determination of the equilibrium position is made by application of the 
principle of minimum potential energy. With equilibrium established, contact stresses 
can be quantified. In order to understand the details of this approach, the principle of 
virtual work and the principle of minimum potential energy must be discussed. 

Given a body in equilibrium, it is desired to describe the response of that body to 
infinitesimal displacements resulting from a system offerees. If each particle in the body 
is described by some generalized coordinates, then the work resulting from these 
infinitesimal displacements is simply the product of the generalized forces acting on each 
particle and the particle's displacement. However, if the particle is in equilibrium, the 
work must be zero since the summation of forces in the x, y, and z directions is zero. 
The infinitesimal displacements and work in this example are referred to as virtual in 
nature. The fact that this work vanishes is referred to as the principle of virtual work. 

The virtual work discussed thus far can be subcategorized as virtual strain energy 
and virtual work done by external forces. From the definition of strain energy, virtual 
strain energy, dU, that results from virtual displacements can be calculated. Since this 
energy is viewed as energy agaifist the bonds between elements, dU is a negative quan- 
tity. The work done by external forces is designated dlV and is simply the summation 
of the product of the external forces and the displacements of the generalized coordi- 
nates. 

Since the principle of virtual work states that the work done as a result of virtual 
displacements is zero. 



dlV-dU = 0 
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Alternately, this can be expressed as, 

n = ()(L-- (D = 0 

where fl represents the system's total potential. 

The above equation illustrates the condition of minimum total potential of a system.. 
This is the foundation of the principle of minimum potential energy. [Ref 12: pp. . 
330-331] This principle states that in a condition of stable equilibrium, the system's 
total potential is stationary. Hence, determination of a system's total potential energy 
and the minimization of that quantity will enable the calculation of the equilibrium po- 
sition. This is the basis for the numerical techniques developed in this study. 

B. CONTACT PROBLEM DESCRIPTION 

The objective of the solution techniques to be developed is to obtain a means of 
calculating the total potential energy and minimize it to determine the equilibrium con- 
dition. To accomplish this, two different models will be created to first solve a simple 
isotropic case, the solution of which is known. In this manner, our models can be vali- 
dated for use in more complicated arrangements. 

Consider contact between a cylinder and an infinite plane as shown in Figure 1. A 
cross section of the contact zone is shown in Figure 2. Let the width of the contact zone 
equal a distance of 2a. .An analytical solution of this problem is available as a result of 
the work done by Hertz. In developing the solution methods, only a very limited region 
adjacent to the contact zone will be e.xamined. The reason for this is that the contact 
phenomenon is a very local one. Figure 3 represents the analytical solution of the 
normal stresses resulting from this contact problem. [Ref 13] As shown, the stresses 
in the foundation diminish very rapidly. The diagram shows that ct, becomes negligible 
at depths less than one half-contact zone (a) and diminishes significantly in less than 
3 half zones. Because of this, it is reasonable to assume that displacements beyond a 
very limited region are negligible in the strain energy calculation of the contact problem. 

,A number of simplifying assumptions will be made for this study: 

1. .As discussed above, displacements beyond a limited region are negligible in strain 
energy calculations. 

2. The foundation is an elastic isotropic material. The cylinder (roller) is rigid. 
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system's degrees-of-freedom and computational time. The penalty method does not re- 
quire additional computational time, but it has been shown to result in less accurate 
solutions in satisfying the contact boundary conditions. Since the contact boundary 
conditions are e.xactly satisfied only when the penalty term goes to infinity, the correct 
choice of the penalty parameter is the key to an acceptable solution. Guerra [Ref. 4] 
supported these claims. 

Because of the limitations of the Lagrange and penalty methods, Bishoff [Ref 5] 
advocated the use of the augmented Lagrange multiplier method to solve finite element 
contact problems. This method is favorable since it avoids the limitations of both 
methods discussed above. The augmented Lagrange multiplier method is similar to the 
penalty method in that the objective function is penalized for constraint violations. 
However, an additional multiplier term is added so the optimum can be achieved by a 
combination of both terms. As stated by Vanderplaats [Ref 6: p. 141], the advantage 
of this method is that the penalty term is not required to grow to infinity to achieve e.xact 
constraint satisfaction. 

The augmented Lagrange multiplier method can be used in numerical programming 
techniques to find the optimum of any functional. Pierre and Lowe [Ref 7] provide a 
detailed analysis of the programming techniques necessarv' in applying this method. 
■Additionally, Vanderplaats [Ref 6: pp. 140-147] provided an excellent discussion on the 
practical usage of this method in numerical techniques. Rothert et al. [Ref S] used a 
numerical programming code based on this method to solve a nonlinear contact prob- 
lem. In this study, an existing numerical optimization routine utilizing the augmented 
Lagrange multiplier method will be used as an integral part of two solution methods 
developed to solve the contact problem. These methods will use different techniques to 
obtain the same functional. The augmented Lagrange multiplier method will then be 
used in a similar fashion to solve each problem. In each case, the optimization routine 
will be used to determine a set of design variables that describes the contact state. 

Following the development and verification of the numerical procedures, this study 
will investigate the response of composite sandwich materials to low velocity impact. 
One of the common failure mode of low velocity impact is delamination. Joshi and Sun 
[Ref 9] studied the impact response of a three layer cross-ply graphite epoxy laminate. 
A correlation was obtained between delamination cracks initiated experimentally and 
maximum shear stress points determined numerically. Sun and Rechak [Ref 10] fol- 
lowed up these findings and found that the introduction of adhesive layers between 
laminae reduced the shear stress distribution thus reducing delamination. Choi, Wang, 
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and Chang [Ref. 11] studied the effects of laminae orientation, ply thickness, and 
stacking sequence on impact damage of graphite epoxy composites. It was determined 
that stacking sequence affects impact damage more than laminae thickness variations. 
.Much of the previous work has focused on the behavior of the graphite epoxy laminate. 
However, there is currently interest in the development of turbine blades constructed of 
sandwich composites. It is therefore beneficial to investigate the response of composite 
sandwich materials to low velocity impact. 
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Figure 1. Roller-foundation assembly 




Figure 2. Contact zone cross section 
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3. Deformations normal to the cross-section are negligible, hence a condition of plane 
strain exists. 

4. The roller is subjected to a vertical distributed load. 

5. The roller-foundation contact is frictionless. 

An important restriction on the minimization problem will be that one body will be 
prohibited from penetrating into the other body. This may seem like an obvious re- 
striction. however a method of mathematically stating this restriction must be discussed. 
Figure 4 shows the deformed and undeformed contact zone. Let () represent the de- 
flection of the roller due to an external force and v(.r^v) represent the vertical deflection 
of the foundation at any point (xj). .At the point of Contact .A, the condition of no 
interference can be expressed as. 



v(O.O) > (5 

■At point B. this condition can be stated as. 

v(r sin 0.0) > (5 — r(l — cos 0) 

The latter condition can be specified at as many points as necessarx' to define this re- 
striction. 

C. NUMERICAL OPTIMIZATION 
1. Optimization Fundamentals 

Before developing the models to be used in this study, one final area must be 
discussed. Once the total potential energy has been calculated, a method of minimizing 
it to find the equilibrium position must be employed. A study of the numerical opti- 
mization technique to be used is required. 

The technique being used in this study relies heavily on the methods of design 
optimization. Design optimization is the utilization of mathematical techniques to 
minimize or maximize a particular value to obtain an optimum solution. The method 
is ideally suited for design. A given design task may have an infinite number of sol- 
utions. However, finding the best solution is a matter of the designer's experience and 
intuition. In the absence of significant experience in a particular field, finding this sol- 
ution may reduce to examining a range of possible solutions by trial and error. Opti- 
mization routines can be utilized to find this solution mathematically. 

Optimization problems can be constrained or unconstrained. For example, it 
may be desired to determine the minimum of the parabola, F{x) = {x — 5)^ + 2. As seen 
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in Figure 5, the minimum is clearly identified at point A. This is an example of an un- 
constrained problem. A constrained counterpart of this problem is: 

•Minimize: F{x) 

Subject to: F(.v) > .5jc -f 4 



As seen in Figure 6, the minimum of the constrained problem is at point B. This is ai 

) 

simple illustration of constrained minimization. The contact problem is a constrained;! 




Figure 4. Deformed contact zone 

The value to be minimized or maximized is referred to as the objective function. 
The parameters to be determined are referred to as design variables. The optimization 
process is an iterative one. The objective function is evaluated. The design vanables 
are varied thus obtaining a new objective function. If the difference is within a certain 
tolerance or meets a certain convergence criteria, the optimum has been obtained. If it 
is out of tolerance, the cycle repeats. 
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Figure 5. Unconstrained minimization 
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Figure 6. Constrained minimization 
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2. Augmented Lagrange Multiplier Method 

The contact problem belongs in the class of constrained problems. There are 
several techniques of solving constrained problems. The technique used in this study 
belongs to a class of solution techniques known as sequential unconstrained minimiza- 
tion techniques (SL’MT). This class of techniques is designed for the general nonlinear 
problem. The fundamental concept behind this approach is that a constrained problem 
is transformed into an unconstrained problem and the objective function is minimized 
using an unconstrained minimization technique. This transformation is accomplished 
by assessing a penalty to the objective function for constraint violations. For example, 
if the design variables are varied in such a way as to enter the region where a constraint 
is violated (i.e., the infeasible region), the objective function would be assessed a penalty 
(i.e., increased). Thus, in order to minimize the objective function, the design variables 
remain within the region where no constraints are violated (i.e., the feasible region). 
[Ref. 6; pp. 121-123] 

There are a number of methods within this class of techniques. They essentially 
differ in the way in which penalties are assessed. The technique used in this study is the 
augmented Lagrange multiplier (.ALM) method. 

Given the constrained inequality optimization problem; 

Minimize; F{X) 

Subject to: < 0, / = 1.2 n 

The .Augmented Lagrangian is defined as. 



m 

A{X, k,p) = F[X) -h + A ] + P\.Si{X) + sff} 

i=i 



where, 

X = vector containing the design variables 
/J., = Lagrange multipliers 
p =penalty parameter 

s, =slack variables which convert inequality constraints to equality con- 
straints 



13 



The first two terms of A(X,/.,p) represent the Lagrangian. f'rom the method of , 
Lagrange multipliers, it is known that minimization of the Lagrangian represents opti- i 
mality. Like simple problems where >. is simply an additional unknown to obtain, in the 
ALM method / is unknown. Hence, a mathematical routine based on constraint values 
is used to select and modify / for successive iterations. [Ref 6; pp. MO- 147] Initial 
selection of this term can have a significant impact of the problem's convergence. 

,\s discussed above, a penalty is assessed to the objective function for constraint 
violations. This feature is apparent by examining the last term. A constraint violation 
results in a positive value for ^(.\') thus resulting in an increase in A. The value of p is 
a scaling term which is sequentially increased throughout optimization. This ensures 
that there is a balance between convergence and numerical conditioning. If p remains 
small, convergence may occur with major constraint violations. If p remains large, 
constraints will be satisfied at the expense of an ill-conditioned problem [Ref 7: p. 
169]. .As with the Lagrange term, selection of this term has a significant impact on the 
outcome of the problem. 

3. Optimizer and One-Dimensional Search 

Thus far. a procedure has been defined which has transformed a constrained 
minimization problem into an unconstrained one. This level of the optimization process 
is referred to as the optimization strategy. The formulation of the modified objective 
function via the augmented Lagrange multiplier method represents a key portion of the 
optimization process. However, there are additional parts of this process that require 
comment. 

With the modified objective function and a procedure for assessing penalties in 
place, a procedure for minimizing the objective function must be defined. This portion 
of the process is carried out by the 'optimizer.' The optimizer's function is to system- 
atically alter the design variables in a manner that reduces the objective function rapidly. 
If A' represents a vector containing the design variables, the following process is used by 
the optimizer to alter X 



A, — A/,_|, -b 



where, 

i s current iteration number 
S = search vector 

k = scalar representing distance traveled in direction S 
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In general, two processes must be accomplished to find the optimum. First, the search 
direction S must be determined by a systematic process. An e.xample of this phase is the 
steepest descent method where the direction of steepest gradient is chosen. Second, the 
scalar k must be determined such that the objective function is minimized as much as 
possible in the search direction of the current iteration. The latter phase is referred to 
as one-dimensional search. [Ref 6; pp. 10-12] 

The optimizer used in this study is the variable metric method. Due to the 
comple.xity of this method, a discussion of their formulation is omitted. A detailed ac- 
count is available in Vanderplaats [Ref 6; pp. 92-93]. The one-dimensional search 
routine used in this study is the golden section method with polynomial interpolation. 
The one-dimensional search portion of the process merely represents a systematic and 
elTicient method of finding the minimum in the chosen search direction. A detailed ac- 
count of this approach is again available in Vanderplaats [Ref 6: pp. 26-49]. 

4. Convergence 

The final point to be discussed relevant to optimization fundamentals is that of 
convergence. Convergence criteria are utilized to identify the optimum solution and 
terminate calculations. There are a number of convergence criteria that can be utilized. 
The most obvious is absolute convergence where the objective functions from two suc- 
cessive iterations are compared. If the difference between the values is within some 
prescribed limit, optimization is terminated. A second method signifying optimality for 
unconstrained problems is calculation of the gradient with respect to the design variable 
vector A'. If this value is approximately zero, optimality has been achieved. This method 
is called the Kuhn-Tucker conditions for unconstrained minimization. Kuhn-Tucker 
conditions are more involved for constrained problems. [Ref 6: pp. 100-101] 

The Automated Design Synthesis (ADS) System used in this study uses both 
these termination criteria as well as relative convergence. Relative convergence is similar 
to absolute convergence except normalized versions of the difference between successive 
iterations is calculated. Again, if a specified tolerance is achieved, optimization is ter- 
minated. 
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III. APPROXIMATE SOLUTION TECHNIQUES 



A. RAYLEIGH-RITZ APPROACH 
1. Background 

The Rayleigh- Rilz method is a method of utilizing the theor\' of minimum po- 
tential energy to solve a given problem. The fundamental concept behind this method 
is that a trial function that represents the deformation field is chosen in terms of un- 
known constants. Second, the system's total potential energy is calculated in terms of 
the trial function. Since the total potential is a minimum at equilibrium, minintization 
enables determination of the unknown constants. The total potential is minimized by 
differentiating with respect to the unknown constants and equating to zero. The result 
is 'n' equations and 'n' unknowns, the trial function constants. [Ref 12: pp. 335-336 ] 

The only requirement of the Rayleigh-Ritz method is that the trial function is 
kinematically admissible. .A kinematically admissible solution is one that satisfies the' 
geometric boundan.' conditions of the system (i.e.. deflection and slope). Other require- 
ments need not be satisfied. .As an example, consider a simply supported beam of length 
L with the origin at the left end of the beam. .A kinematically admissible solution to 
describe the beam's one dimensional deformation from vertical loading in the y direction 
is, 

>•(.0 = sin( ) 

where a„ represents the coefficient to be determined. 

Deflection boundary conditions at x = 0 and x — L have been satisfied. Obviously, an 
increased number of terms in the trial function will yield a far more accurate solution. 
A Fourier sine series would be a reasonable selection in this case. 

.Although the Rayleigh-Ritz method does not stipulate numerous requirements 
on the trial function, sensible choices of trial functions will increase solution accuracy 
significantly. For example, consider the same simply supported beam with the origin 
now at the center. A kinematically admissible function is, 

j'(jt) = a„ sin( ). 
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Naturally however, due to the placement of the origin in this problem, an even function 
is a much for sensible selection for a trial function. A more appropriate selection would 
be, 

>-(.v) = a„ cos( ^ ). 

The latter point concerning sensible selection of the trial function will be discussed in 
detail throughout this study. 

2. Application of the Rayleigh-Ritz Method to the Contact Problem. 

A brief overview of the method to be developed is in order. .Application of the 
Rayleigh-Ritz method necessitates the selection of the appropriate trial function in terms 
of unknown coelTicients. A discussion of the physical nature of this problem as well as 
the desired features of the trial solution is required. 

The theorv’ of elasticity relationships will be applied using the trial function to 
obtain the system's strain energy in terms of unknown coefficients. The system's total 
potential energy will then be minimized utilizing the optimization techniques discussed 
in Chapter II. The design variables are the unknown trial function coefficients. With 
the coelTicients determined, the displacement is known for all points enabling the stress 
to be determined throughout the body. Figure 7 is a flow chart of the procedure to be 
utilized. The post-processing procedure shown is simply the calculation of the stresses 
using the now determined coefficients. 

3. Trial Function Selection 

In order to choose an appropriate trial function, it is necessarv' to have an 
understanding of the physical phenomena to be modeled. Consider the roller-foundation 
system shown in Figure 1. Two trial functions are needed to model horizontal and ver- 
tical deformation. There are a number of important characteristics that should be in- 
herent within the trial functions. These are outlined below: 

1. As seen in Figure 3, <j, and are equal and compressive at the point of contact. 
Additionally, <7, decays more rapidly than as distance from the contact point in- 
creases. 

2. Taking the origin at the point of contact as shown in Figure 2, the u-deformation 
takes the form of an odd function. 

3. The v-deformation takes the form of an even function (i.e. symmetric about the 
origin and greatest at the origin). 

4. To satisfy the Rayleigh-Ritz method requirements, the boundary conditions must 
be satisfied. In this case, this requires u and v deformation to be zero at the far 
boundaries. 
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Figure 7. Rayleigh-Ritz method applied to contact problem 
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5. The characteristics outlined in item 4 satisfy the requirements of this method. 
However, since this study will be a stress analysis, it is also desired that the stresses 
also reflect the physical phenomena. Since a, and are functions of e. and c, 
and should also exhibit certain characteristics. Referring to Figure 8 for dimen- 
sions and the coordinate system, it is desired that r, and r, equal zero at x = L 2 and 
y= H. This will ensure that stresses are zero at the boundaries. To satisfy the re- 
quirements of item 1 above, c, and should be equal and negative at the point of 
contact and decrease in magnitude as the distance from the point of contact in- 
creases. 




Figure 8. Contact zone 

With the above guidelines in mind, the trial function can be selected. The trial 
functions chosen for this study are composed of a series of terms of the general form: 



u(x^) = aJH-yr(j--xf (3.1) 

sixy) = b„{H-y)\\-xf (3.2) 

These expressions were carefully chosen and represent a compromise due to the diffi- 
culties of satisfying all boundary conditions with the physical phenomena of this prob- 
lem. 
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From the above expressions, it is immediately obvious that the geometric 
boundary conditions at x=L 2 and y=H have been satisfied. This satisfies the re- 
quirements of the Rayleigh-Ritz method. In addition, there are a number of important 
characteristics that illustrate the advantage of this selection; 

1. £, and £j, are negative thus simulating a compressive environment in the vicinity of 
the point of contact. The importance of this is obvious. If normal strains were not 
negative, the resulting requirement would be for the coefficient to be less than zero 
to simulate compression. It is obvious that this would result in deformations op-^ 
posite to that which was desired by the choice of the trial function. 

2. This selection for deformation fields has the important characteristic of decreasing; 
deformation as we move away from the point of contact. .Also note that defor- 
mation is maximum at the point of contact. 

3. The exponents a.b.c. and d can be varied to simulate subsurface stress fields. For' 
example, if the analytical solution indicates a large y gradient for <r, at x = 0, the' 
objective would be to increase the rate at which e, decreases as the distance from 
the point of contact increases. This would be easily simulated by raising the value' 
of a. If this change had a detrimental effect on the behavior of the exponents- 
of the vertical deformation could be varied to restore the solution. 

It is important to note that this selection is not without compromise. The most 
significant limitation of the trial functions is with regard to the horizontal deformation 
u. Physically, it is expected the behave as an odd function as discussed above. 

However, in this selection of trial function, a positive value of ir(.rj-) exists at the origin. 
This is contrary' to the physical behavior of the problem and will lead to some error. 
However, considering that the magnitude of this deformation in the elastic range is 
small, this error is believed to be limited. .Another consequence of this compromise is 
the existence of non-zero shear strain at jc = 0. 

.Another less severe limitation is a restriction on the order of the exponents in 
order to maintain zero stress at the boundaries. Since ct, and are combinations of £, 
and £^, both normal strains must be zero at the boundaries to ensure that stresses are 
zero at these locations. Since the normal strains are first derivatives, this requires that 
b and c are at least equal to 2. 

It is worthwhile to note that most of the considerations discussed above far ex- 
ceed the requirements stipulated by the Rayleigh-Ritz method. The objective has been 
to utilize trial functions that closely match the physical nature of the problem in an ef- 
fort to maximize accuracy. 

The specific trial functions used in this study are listed below. The horizontal 
deformation was assumed to be; 
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(3.3) 



n=\ 

where // and L represent the height and length of the bearing foundation, respectively. 
Summation was done for n equal 1 and 4. The vertical deformation was assumed to be; 

m 

v(.rj-) = y^b,{H -y)-{ (3.4) 

A2=l 



As discussed previously, manipulation of the exponents enables the trial function results 
to he matched with the analytical solution. As will be illustrated in the results, the ex- 
ponents chosen in the above functions achieve this goal sufTiciently. 

With the deformations chosen, the stresses and strains can be calculated for use 
in the calculation of the foundation strain energy. These values are shown below: 



m 

,, = y-2a,{H-yf^^\j--x) 

n=\ 


(3.5) 


i=\ 


(3.6) 


(l+v)(l-2v) 


(3.7rt) 


(1 +v)(l -2 v) 0£y+V£j 


(3.76) 



where E and v are Young's modulus and Poisson's ratio, respectively. For the shear 
stress and strain, 



iH- = ^ _ (4 -f n)a,{H-yf^''\ j- - xf (3.8) 

n=] 

m 

= ^ - ( 1 + n)bJ,H - y)\ j- - xf (3.9) 

(=1 
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(3.10) 




du cv 




(3.11)1 



where G is the shear modulus. 

Using the above quantities, the strain energy L' can be calculated. From the definition, 
of strain energy applied to a two dimension deformation field, 



Because of symmetry about the origin, strain energy can be calculated for half of the 
domain and doubled. With strain energy calculated, the total potential for the system 
can be found from. 



where, 

F = e.xternal force per unit length applied to the roller 
(5 = vertical distance traveled by the roller. 

The quantity F3 represents the work done by the roller on the bearing foundation. 

B. FINITE ELEMENT APPROACH 
I. Total Potential Derivation 

The finite element method can be employed to solve the contact problem. A 
finite element mesh can be constructed to appro.ximate the behavior of an elastic foun- 
dation subjected to line contact loading from an rigid roller. The resultant interaction 
between the foundation and roller enables the calculation of the foundation's strain en- 
ergy and the subsequent calculation of the total potential energy. By again utilizing the 
optimization techniques discussed in Chapter II, the equilibrium position can be deter- 
mined. Thus the contact stresses can be calculated throughout the body. 

The objective is to derive a means of calculating the total potential of the system 
by application of the finite element technique. Total potential energy is defined as, 




(3.12) 



n = u - f 



(3.13) 
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n = U-F3 



( 3 . 14 ) 



where, 

i' = iniernal strain energy 

F= external force per unit length applied to the roller 
3 = vertical distance traveled by the roller. 

The strain energy of the system can be found from. 



L = 



-) (^,x4"x ^ xy 'f xy)^^^ 



a 



( 3 . 15 ) 



where Q represents the problem domain. This can be expressed in matrix form as. 



where. 



i' = 






Q 



{e}^= {£, Cy y^y} 
{^X <^y ^Xy} 



( 3 . 16 ) 



On the element level, 



t 



U = 









( 3 . 17 ) 



where i represents analysis of the element. 

The stress matrix can be expressed as, 

{^} = [/)]{£} ( 3 . 18 ) 

where [T>] represents the material property matrix. 
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For a condition of plane strain, the stiffness matrix can be expressed as, 



[Z)] = 



F 

( 1 + V’)( 1 



2v) 



I - V V 0 

V 1 — V 0 



For a plane stress condition. 



[^] 




1 

V 

0 



V 0 

1 0 



(3.19a) 



(}.\9h) 



The development of this technique will use linear triangular elements. The method, 
however, can be applied to any type element. For linear triangular element, the defor- 
mations take the following form: 


u = 


(3.20) 


V = //]Vj -(- ^ 2 ^ 2 ^ 3^3 


(3.21) 


where the shape functions H, are defined as: 




- [(-'^iV’3 - + Cv ‘2 >’ 3 )-'^ + (•'^3 ■) 


(3.22) 


- •^l>'3) + 0^3 - + (-'^1 - -'^3 )t] 


(3.23) 


^3 - [(•^ 1>2 - -'■ifl) + CV'I ->• 2 )-'^ + (-’^2 - -^i)t] 2J 


(3.24) 



where, 

= coordinates for node / 

A = element area [Ref 14 ]. 

The strain matrix can be expressed as, 
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{£} = 



d 

cx 

0 



d 

cv 



d d 
cv cx 



Substituting equations (3.20) and (3.21) into equation (3.25) yields, 



{£} = 



cH 



cx 

0 

cH, 



0 

c//, 

cv 

cH^ 



clL 



cx 



0 



0 

cH, 



cH 



3 



CV 

cH-, c//i 



O’ 



cx C V 



cx 



CX 

0 

C//3 

cy 



cH 



3 



O' 

cH 



3 



cx 



U2 

V’2 

"3 



In abbreviated notation, equation (3.26) is expressed as, 

{£} = [5]{^ 

For the linear triangular element, [fi] reduces to, 






>2 ~ >'3 0 }'3 ~ d'l >’l ~ y'l ^ 

0 -Vj - .V2 0 -V, - .V3 0 ,V2 - -V, 

-'^3 -'^2 >’2-d’3 •^l--^3 J’ 3 - 3’1 -'’2--^l >’l “>'2 



(3.25) 



13.26) 



(3.27) 



;3.2S) 



Returning to the element strain energy calculation, equation (3.17), 



U = 



\{z]\o]dQ. 



and substituting (3.18) into (3.17), 



U = 



‘'q' 



(3.29) 
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Substituting equation (3.27) yields, 



U = 






"O’ 



r = ^{./}^[5]^[Z)][5].q/{^ (3.30) 

where i = unit depth. Defining the element stiffness matri.x , 

[A-] = [5]^[f)][5].-l/ (3.31) 

then the strain energy per element equals, 

r = ^{./}[A’]{./} (3.32) 

The element stiffness matri.x can be expanded into the global stiffness matrix. 
With strain energy now determined, total potential can be determined from equation 
(3.14). 

2. Optimization and Static Condensation 

.A.S with the previously developed model, the augmented Lagrange multiplier: 
method will be utilized to determine the equilibrium condition via the theorem of mini- 
mum potential. The objective function is again the total potential. In this case, how- 
ever, the design variables are the nodal deformations, u, and v,, for non-fixed nodes. 
Constraint equations are developed in a similar manner as discussed in Chapter 11 to: 
ensure that one body does not violate space occupied by the other. 

Since the nodal deformations are represented as the optimization design vari- 
ables, the number of design variables will equal twice the number of non-fixed nodes. 
For a simple mesh, a direct application of this procedure will likely yield accurate results. 
However, it is known that the accuracy of the optimization routine declines as the: 
number of design variables increases. Hence, for complicated finite element meshes, 
solution accuracy will be adversely affected by the large number of design variables. 
Therefore a procedure must be adopted to eliminate the need for assigning design vari- 
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ables to nodes where information is not necessar\' for evaluating a solution. This pro- 
cedure is known as static condensation. 

Static condensation has been utilized by References 2 and 3 in an elTort to re- 
duce computer computational time. The idea behind static condensation is the reor- 
ganization of the global stilTness matri.x. A finite element problem can be expressed as, 

= {F} (3.33) 



where. 

[^;] = the stilTness matrix 
[u] - the deformation vector 
{F} = the force vector 



It is desired to reorganize this system of equations into the following: 

^11 ^12 
^21 ^22 




(3.34) 



The vector n, contains the essential nodes while vector n, contains non-essential nodes. 
Essential nodes are those where boundary’ conditions are applied and nodes that are as- 
signed optimization design variables. By matri.x manipulation, 

{w2} = -[A'22]-'[A.'3,](n,}. (3.35) 

Therefore, the displacement vector can be expressed as, 

{“i) 

-[A'22]‘'[A2,]{u,} 

where 1 represents the identity matrix. 




[/] 

-[^22r'C^'2,] 



U‘\ 



(3.36) 



Substituting this equation into the global counterpart of equation (3.32), 




[/] 

[A^22r’C^2l] 






[/] 

[^22]"'[^'2l] 






By defining the reduced stiffness matrix [A], the above reduces to. 



(3.37) 
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(3.38) 



= y («,}''[ 

As discussed at the beginning of this section, the ADS design variables are as- 
signed as the horizontal and vertical deformations at all non-fixed nodes. With the in-- 
tegration of static condensation, design variable assignments are further restricted tO) 
non-lixed. non-condensed nodes. The procedure is now in place for calculation of strain i 
energy and total potential energy. In summary, a flow chart of the solution procedure' 
utilized in this chapter is shown in Figure 9. 
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Figure 9. Finite element method applied to contact problem 
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IV. RESULTS AND DISCUSSION 



A. PROCEDURE VALIDATIONS 
1. Rayleigh- Ritz Method Results 

In Chapter III an approximation technique was developed to solve the contacti 
problem via the Rayleigh-Ritz method. .As discussed, two trial functions that approxi-- 
mated the deformation field were selected in terms of unknown coelTicients. The hori- 
zontal deformation was assumed to be; 



Using the above deformation fields, theory of elasticity relationships and the definition 
of strain energy were employed to obtain an expression for the total potential energy of 
the system shown in Figure I. Numerical minimization techniques were then employed 
to determine the equilibrium condition and the contact stresses. 

To illustrate the application of this method, an isotropic material with the fol- 
lowing properties was selected; 



.As stated in Chapter II, this problem was selected for development of this technique 
because an analytical solution is available as a result of the work done be Hertz. It is 
desired to use this analytical solution to choose a roller size, load, and problem domain 
that can be used to accurately simulates the contact phenomenon. 

Contact stresses as well as the size of the surrounding region of influence are 
strongly affected by the size of the contact zone. Naturally, as the size of the contact 



m 






The vertical deformation was assumed to be; 



m 




E = 200 GPa 
v = 0.3 

G = 76.9 GPa 
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zone increases, the load is distributed over a larger area and contact stresses decrease. 
The e.xtent of the affected subsurface zone also decreases. From the analytical solution, 
the contact zone size is defined by the e.xternally applied force, the material properties 
and the diameter of the roller [ Ref 13 ]. Hence for a given material, the roller size and 
the e.xternal load define the contact zone size and the resulting stresses. 

Figure 2 defined the width of the contact zone as 2a. Figure 3 shows the ana- 
lytical solution of the contact problem. This figure shows the decrease of the normal 
stresses as a function of half-contact zones (a) awa> from the contact point. As shown. 

decreases more gradually than o. . Therefore the decay of is the limiting factor in 
defining the domain beyond which strain energy contributions are negligible. From 
Figure 3, it is estimated that the contact phenomenon can be accurately modeled by 
exantining a region equal to approximately five half-zones (5a). 

Since a numerical integration technique was used to perform the double inte- 
gration required by Equation 3.12, the dimensions were selected for numerical conven- 
ience. Referring to Figure 8, height H and length L were selected as 1 and 2 meters, 
respectively. Due to the problem's symmetry, half the foundation was analyzed. This 
enabled the double integration to be conducted between the limits ofO and 1. Since the 
foundation height H has been set to 1 m, it is desired to have this distance equal to 5 
contact zones (5a) as described above. 

Using the analytical solution, a load and roller diameter were selected that cre- 
ated a contact zone such that the foundation height H was equal to a distance of 5a. 
The only additional restriction was that the resulting contact stresses remained within 
the elastic range of the material. Yield stress was assumed to be 300 MPa. The load 
and roller radius combination used in this study are; 

Load: 90 .MN 
Radius: 75 m 

The values were obtained using the analytical solutions found in Reference 13. The 
latter dimension is unrealistic. However, as stated above, it is a result of the selection 
of base dimensions in the interest of numerical convenience. This model simply repres- 
ents a scaled-up examination of the base material in a small region adjacent to the con- 
tact zone. 

Comparisons of the approximated contact stresses with the analytical solution 
along the axis of symmetry are shown in Figures 10 and 1 1. Both figures are normalized 
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graphs of the stress distribution. As shown in Figure 3, the maximum stress occurs at 
the point of contact. As shown, this stress is equal to a, and at the point of contact. 
Figure 10 is a comparison for . This figure shows a stress distribution that closely 
approximates the analytical solution. Figure 1 1 shows the analytical and approximate 
stress distribution for a,. It is apparent that this method over approximates this stress. 
.\s stated in Chapter 111, one of the benefits of this trial function is the ability to change 
the exponents to match analytical solutions with approximate solutions. A brief expla- 
nation of the choice of exponents used in this solution technique follows. 

While selecting exponents, it must he understood that ct, and are both func- 
tions of horizontal and vertical deformations. Therefore changing the exponent of one 
deformation to alter the stress in one direction will influence the behavior of the other. 
In the case of Figure 11. it appears that a modification is required. ,\n increase in the 
exponent of the y-portion of the horizontal deformation function seems appropriate to 
increase the rate at which 0 , decays. However, this action will have the undesirable ef- 
fect of decreasing the rate of decay of It is possible to counter by decreasing the; 
exponents of the vertical deformation. However, as stated in Chapter HI. in order to 
maintain a zero stress boundary condition at x=L 2 and y=H the exponents of all 
terms must be greater than or equal to 2. Thus a compromise must be reached. It is 
believed that since is the dominant term, priority should be placed on ensuring is 
as close as possible to the analytical solution. .Accordingly, a decision was made to ac- 
cept the over estimation of cr, as shown in Figure 11. In this case, the estimation of a,, 
is conservative. 

The contours of the Rayleigh- Ritz solution are shown in Figures 12 through 15.. 
It has been determined that this method can accurately predict contact stresses resulting: 
from line contact between a roller and flat plane. Table 1 compares the approximated! 
contact stresses and the analytical results at the point of contact. 



Table 1. RAYLEIGH-RITZ RESULTS AT CONTACT POINT 





Current .Model 


Analytical Solution 


c7. (.MPa) 


285.6 


289.7 


a, (MPa) 


301.9 


289.7 
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Figure 1 1. Comparison of <r, in a roller contact problem 
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Figure 12. Stress contour of from Rayleigh-Ritz method 
2. Finite Element Method Results 

In Chapter III, an approximation method was developed to solve the contact 
problem using the finite element method. A method of calculating a system's strain en- 
ergy and the total potential energy was investigated. In addition, the use of static 
condensation to improve optimization efficiency was described. As discussed, the nu- 
merical minimization techniques were again utilized to determine the equilibrium posi- 
tion. A means of employing these techniques to evaluate the contact phenomenon was 
introduced. In this section, a simple contact problem will be investigated to validate the 
algorithms used to calculate strain energy and those used to implement static 
condensation. With confidence in these algorithms, the more complex roller-foundation 
problem will be examined. 

a. Two Thin Plates in Contact 

The procedure developed was first validated on a simple contact problem 
the solution of which was known. In this example, two thin bodies in plane stress were 
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Figure 13. Stress contour of from Rayleigh-Ritz method 

studied. As shown in Figure 16, one body, restrained on one edge and subjected to a 
horizontal load, comes in contact with a second body rigidly supported on three sides. 
The finite element model developed to solve this problem is composed of 14 linear tri- 
angular elements as shown in Figure 17. The objective of the fortran program developed 
to solve this problem was to calculate the total potential energy of the system using the 
finite element technique and the method of static condensation. With this accomplished, 
the equilibrium position can then be found via the augmented Lagrange multiplier 
method. 

The objective function for this problem is the total potential energy. 
Equation 3.14. Referring to Figure 16 and 17 the constraints imposed on the system are 
expressed as: 



Wg < 0.005 + Wn 
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Figure 14. Strain contour of from Rayleigh-Ritz method 

u<) < 0.005 + W|2 

where n represents the horizontal deformation of node i. 

.As shown in Figure 17, seventeen nodes were used to model the system. 
Each node has been assigned horizontal and vertical deformation variables. .Accord- 
ingly, the degree of freedom for this system is twice the number of nodes. .As discussed 
in Chapter III, static condensation requires the identification of essential nodes and 
non-essential nodes. Essential nodes are those nodes where .ADS design variables are 
assigned and boundar\' conditions are applied. Referring to Figure 17, node 3 is the 
point of load application and must be assigned a design variable. Nodes S, 9, 11, and 
12 are assigned design variables in order to define the constraint equations described 
above. .After eliminating all fi.xed nodes from consideration, nodes 2, 5, and 6 are the 
only candidates for condensation. .As discussed in Chapter III, the global stiffness ma- 
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Figure 15. Strain Contour from Rayleigh- Ritz method 

trix is rearranged according to Equation 3.34. For this problem, the vector containing- 
the non-essential nodal information, w,, is arranged as follows; 

{ W 2)^={«2 1’2 ^^5 I's »6 ^ 6 } 

where u, and v, represent the horizontal and vertical deformation of node /, respectively. 
Strain energy and total potential energy were calculated according to Equations 3.38 and 
3.14. The latter was minimized using the augmented Lagrange multiplier method de- 
scribed in Chapter II. 

The solution obtained from this simple problem were compared with the 
results obtained from Y.W. Kwon and J.E. Akin [ Ref 15 ] and are shown in Table 2. 
The solutions were in agreement. It was concluded that a satisfactory procedure was in 
place to solve the more complex roller-foundation problem. 
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Figure 16. T«o thin plates in contact 



Table 2. FINITE ELEMENT RESULTS: TVVO PLATES IN CONTACT 



Load (N) 


Deformation 


Current .Model 


Reference 15 


l.OxlCH 




No Contact 


No Contact 


1.0,\10» 


u, 


.503x10-2 


.503x10-2 




«I2 


.336x10"^ 


.319x10"' 


1.0x10’ 


“9 


.702x10-2 


.734x10-2 




“12 


.201x10"^ 


.235x10"^ 



b. Roller- Foundation Contact Problem 

A finite element grid composed of 512 linear triangular elements was con- 
structed to model the roller-foundation assembly shown in Figure 1. Because of the 
symmetr>‘ of the problem, one-half of the foundation was modeled. A refined mesh was 
constructed in the vicinity of the point of contact. The mesh is shown in Figure 18 
where the origin represents the point of contact. The domain dimensions are similar to 
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Figure 17. Finite element mesh for t«o plates in contact 
those chosen in the Rayleigh-Ritz method discussed in Part A of this chapter. The roller 
radius was chosen as 75 meters and the half-domain dimensions are 1.40 x 1.40 meters. 
As discussed in Part A, these dimensions represent an analysis of the region immediately 
adjacent to the contact zone and are a result of the local nature of the contact problem. 

Constraint equations were constructed according to the discussion of 
Chapter II Part B. Boundary conditions were imposed in the follouang manner: 

• Horizontal and vertical deformations were prohibited on the remote mesh bound- ' 

aries (i.e., u(1.40,y’) = i;(l .40,>') = 0 and u(jt, 1.40) = v’(x,1.40) = 0). ’ 

• Horizontal deformation was prohibited on the axis of symmetry' (i.e., u(0,>’) = 0). 
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Figure 18. Finite element mesh for roller-foundation problem 
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Since the system has 2S9 nodes, the resulting 578 degrees of freedom ne- 
cessitated the utilization of static condensation. Referring to the contact surface, the 
deformation variables that correspond to nodes on this surface are required for con- • 
straint equations. The nodes that comprise the other three borders of the domain are 

il 

all subject to boundary conditions. Consequently, the interior nodes of the mesh are the ' 
nodes that are candidates for static condensation. In this model, all interior nodes were. 

condensed. The original system was reduced from 578 to 128 degrees of freedom. After * 

! 

application of the boundary conditions, there were 46 possible deformations, a suffi- 1 
ciently small number of design variables for the optimization algorithm. One additional | 
design variable was used to represent the distance of travel by the roller. This value is i 
needed to calculate the work done by the roller on the bearing foundation. j 

As before. Equation 3.38 was used to calculate the system's strain energy. 
Following calculation and the subsequent minimization of the total potential energy, a 
post-processing procedure was followed to determine the contact stresses. The output . 
of the optimization routine represents the nonzero components of the {t/,} vector. In * 
order to calculate stresses throughout the body, the remaining deformations contained’! 
within the condensed vector {;/,) must be determined. This vector is calculated directly I 

i 

using Equation 3.35. With deformations known throughout the domain, strains can be ) 
determined by applying Equation 3.27 to each element. The subsequent application of ! 
Equation 3.18 enables determination of stress for each element. 

To illustrate the capability of this method, an isotropic material with the 
following properties was chosen: 



E = 240 GPa 
V = 0.3 
G = 92.3 GPa 
Load = 90.0 MPa 

Figure 19 shows the deformation resulting from the loading. For clarity, the defor- 
mations have been magnified 100 times their original values. Comparisons of the stress 
distributions with the analytical solution along the axis of symmetry are shown in Fig- 
ures 20 and 21. These figures are similar to Figures 10 and 1 1 and represent normalized 
versions of contact stresses. As shown in these figures, this method is a good approxi- 
mation of the stress distribution in the foundation of a loaded roller bearing. If the mesh 
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was more refined near the contact zone and the domain extended further, the agreement 
between the numerical and analytical solutions would be better. Figures 22 and 23 rep- 
resent normal stress contours of this problem. Figures 24 and 25 show normal strains. 
Table 3 shows a comparison of the results of this model and the analytical solution at 
a selected element in the region of contact. 



Table 3. COMPARISON OF STRESSES NEAR THE POINT OF CONTACT 





FE.M Solution 


.Analytical Sol- 
ution 


a, (MPa) at x = 0.0137. y = 0.0273 


281. 3 


272.0 


CT, (MPa) at \ = 0.0137. y = 0.0273 


315.7 


313.8 



B. APPLICATION 

The preceding section illustrated that the contact problem can be accurately simu- 
lated using the methods developed in Chapter HI. It is the objecti\ e of this section to 
show how this approach can be applied to a contact problem in a composite plate sub- 
jected to low- velocity impact. 

A multi-ply laminate model has been constructed to investigate the response of 
composite materials to low velocity impact. It has been found that composite bodies 
subject to impact damage commonly fail due to delamination. Sandwich composites are 
currently being considered for use as turbine blades. It would be beneficial to acquire 
an understanding of the behavior of sandwich materials to impact damage. 

In order to accomplish this task, a clamped composite beam similar to the one de- 
picted in Figure 26 has been modeled. The beam length is 25 cm. Beam thickness is 2.5 
cm. Because of symmetry, half the beam was modeled with 256 bilinear elements. As 
shown in Figure 27, the mesh is refined near the point of contact, the origin of the mesh. 

The major assumptions of this model are that the beam is in a condition of plane 
strain and that the dynamic effect of the impact can be neglected. Reference 10 ap- 
proximated the loading resulting from low velocity impact as a sinusoid. In this study, 
the peak load will be examined to study the maximum normal and shear stresses. .Ac- 
cordingly, the stress distributions obtained from this study will represent a 'snap-shot' 
in time of the response of the body at the instant of maximum loading. 

The finite element program developed to solve this problem is sufficiently flexible to 
alter the material stiffness matrix [Z)] shown in Equation 3.18 during construction of the 
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Figure 20. .Analytical solution vs. appro.ximate from finite element results 
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Figure 22. Stress contour from finite element results (increment 0.01 GPa) 
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Figure 23. Stress contour a, from finite element results (increment 0.02 GPa) 
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Figure 24. Strain contour from Finite element results (increment 0.000035) 



49 




Figure 25. Strain contour e, from finite element results (increment 0.00001) 
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Figure 26. Clamped composite beam 
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Figure 27. Finite element mesh for clamped beam model 

finite element global siilTness matrix. Therefore by defining the layup for the laminate, 
the lamina stiffness matrices can be varied from element to element to accurately model 
the behavior of the body. This enables a variety of laminate layups and lamina thick- 
nesses to be studied. 

The sandwich materials used in this study are composed of an isotropic interior 
material and orthotropic exterior laminae. Since a condition of plane strain was as- 
sumed, the isotropic material stiffness matrix is given by Equation 3.19a. For the 
orthotropic e.xterior laminae, the material stiffness matrix is given by Equation 4.1 
[Ref 16]. 
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* ~ ''*23 *^32 ^’l2 '^32'^I3 

^' 2 ^ 3 ^ £|£3 7 - 

f|2 + V32V’|3 1 ~ ^’|3''’31 

E,E,T E,E,T 

0 0 G,2 



(4.1) 



where. 

1 - - V:3V3: ~ ^’si^n ~ -V’:|V3;V|3 

£,£,£3 

£, = Young's Modulus in direction 

v ,= Poisson's ratio for lateral contraction in _/'" direction resulting from loading in 
the E direction 



The derivation of the total potential energy calculation in Chapter 111, Part B was 
done using linear triangular elements. Since bilinear elements were used in this model, 
calculation of the element stiffness matrix was more computationally intensive. Indi- 
\ idual entries of the element stiffness matrix were obtained from Reference 17. Other- 
wise, the calculation of total potential energy was identical to the procedure outlined in 
Chapter III. 

There were two groups of boundarx' conditions applied to the problem. Along the 
clamped edge, deformation was prohibited. In addition, horizontal deformation was 
prohibited along the axis of symmetrx'. As before, the application of static condensation 
requires the identification of essential nodes the information of which is contained within 
the {U|} vector. In addition to the essential nodes associated with the above boundarx 
conditions, the nodes along the contact surface are needed for the constraint equations. 
.All other nodes were condensed out. 

To illustrate this application of problem solving, sandwich material with the fol- 
lowing properties were studied: 

Exterior Laminae 
£„ = 170 GPa 
E^= 11.8 GPa 
G ,3 = 5.2 GPa 
v,2 = 0.33 
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Isotropic Core 
E = 2.24 GPa 
G = 0.84 GPa 
V = 0.35 

The beam was loaded by contact with a 10 cm radius ball. The transmitted force . 
was 250 N. Some une.xpected trends were observed in the equilibrium position deter- 
mined by the optimizer. By examining the deformations along the axis of symmetry', a 
gradually decreasing trend in deformation moving away from the point of contact was' 
interrupted in lamina of significantly decreased stiffness. It is believed that this difficulty 
resulted from an inability of the optimization routine to approximate deformations 
through regions containing very dilTerent orders of strain energy. In spite of these dif- 
ficulties, some critical information was obtained from this program. As discussed in the' 
beginning of this study, one of the greatest dilTiculties of the contact problem is the de- 
termination of the size of the contact zone. Fortunately, the size of the contact zone can 
be readily determined by examining the output from the constraint equations. By com- 
paring the ball radius (r) and the distance between the ball center and the node (r'), it 
can be determined if a node is in contact. This condition is illustrated in Figure 28. The 
distance r' to the i„ node is given by the equation: 

r' = \ (i" ~ d -I- + xj 

If r' is greater than r, the node is not in contact with the body. 

With the extent of the contact zone known, the solution to this problem was ob- 
tained by applying the contact boundar>’ condition to a direct finite element program. 
Since the validity of the optimization program was in question, this method was applied 
using a 0-90-0 layup similar to one used in a study conducted by Sun and Rechak [Ref 
10]. The solution obtained from the current approach was very close to the other re- 
sults. 

With confidence in the procedure, it was desired to examine the behavior of this 
model to various layups. The objective was to illustrate how this solution technique can 
be used for meaningful research. The study conducted by Sun and Rechak analyzed 
methods of reducing the likelihood of composite failure due to delamination. Of par- 
ticular concern is the magnitude of shear stress distribution and tensile stress in the y- 
direction, the two predominant causes of delamination failure. Using materials with the 
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thickness of the exterior layer increases, the maximum value of this stress increases. 
However, this stress is always compressive at the interface with the core. By comparing 
the magnitudes of the shear and transverse normal stresses at the interface, it would 
appear that if delamination was to occur at this interface, it is more likely to be caused 
by high shear stresses. 

Comparisons of the bending stresses at cross section A are shown in Figures 39, 40, 
and 41. These figures show that as the thickness of the exterior layer increases, a non- 
linear stress distribution intensifies in the layer closest to the contact surface. This trend 
would indicate that beam theory is unsuitable for estimating bending stress through this 
lamina. This nonlinear behavior is local to the contact zone. Figure 42 shows the 
counterpart for Figure 41 at cross section C. The stress distribution in the contact layer 
is approximately linear. Another significant observation can be made by e.xamining the 
bending stress graphs. As the thickness of the exterior layer opposite to the contact 
layer increases, the stress distribution within this layer transforms from purely tensile 
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behavior to compressive-tensile behavior. This is significant because a bending crack 
initiated by tensile stresses tends to propagate to the core interface and cause delami- 
nation. The presence of compressive stress within this layer will tend to slow the growth 
of this crack toward the core. 

Thus far symmetric layups have been studied. To analyze how beams with asym- 
metric layups respond to contact loading, two beams with the following designations 
were studied: 0(3)-ISO(6)-0(7) and 0(8)-ISO(6)-0(2). The first designated layer is the 
lamina closest to the contact surface. Figures 43 and 44 show the shear stress distrib- 
utions for these tw’o layups. As before, the ma.ximum shear stress is relatively unaffected 
by the different layups. As was the case for the symmetric beams, a thicker exterior layer 
close to the contact zone results in a more gradual transition of shear stress into the 
core. The result is a lower shear stress at the interface for the 0(8)-ISO(6)-0(2) case. 
.-\s seen in Figures 45 and 46, the trends for the transverse normal stress, a^, were the 
same as those found in the symmetric beams. Deflection was less for the 
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properties outlined above, a sandwich composite beam with outer fibers aligned to 0 
degrees was first studied. Since the finite element model is composed of 16 layers, the 
beams studied will be described with the number of finite element layers in parenthesis 
following the layer description. For example, a 0(3)-ISO( 10)-0(3) beam is composed of 
10 isotropic core layers within 3 layers of material with the fibers oriented at 0 degrees 
on the top and bottom of the beam. 

Three symmetric layups of varying core thickness were initially studied. The beams 
have the following designations: 0(3)-1SO(10)-0(3), 0(4)-ISO(S)-0(4), and 

0(5)-ISO(6)-0(5). Figure 29 shows the deformed 0(3)-ISO( l0)-0{3) beam with defor- 
mations magnified 100 times. Stress contours for this beam are shown in Figures 30, 
31, and 32. Figure 30 shows the shear stress contour for the loaded condition. This 
stress is of particular concern since delamination, a common failure mode for compos- 
ites, is commonly initiated by high shear stresses or tensile transverse normal stresses. 
.As the figure shows, a very high stress gradient is present near the contact zone. .As the 
distance along the beam increases away from the contact zone, the magnitude of the 
gradient decreases until the cross sectional shear stress distribution becomes parabolic. 
The transverse normal stress is also concentrated around the contact zone. 

Figures 33, 34, and 35 show the cross sectional shear stress distributions for the 
three symmetric beams described above. Three separate cross sections are shown on 
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C,ch .raph .he locauons or which are indicated in Figure 27. The vertical dashed lines 
on each erapl, identifv the lamtna tnterraces. These ftgures show that all stgn.hcant ac- 
::t :onLed to the lanuna closes, to the contact zone. 1. ,s also e. .. tha 

,s more gradual at the inierfatp with the core. The result is a reduction in 
stress at the interface for layups with thicker exterior lamina. It is also noteuor , 
rhe maximum shear stress occurs at cross section B vice cross section A. 

Figures 36 37 and 38 are graphs for the transverse normal stress. , , 
jriel llyups. These graphs show an increase in . as - .hic^ness of ^ 
■t or lavers increases. The increased thickness of these layers produces a strong 
hence beam dehection and contact zone size are reduced. ’ 

stresses increase. At cross section C. some tensile transverse stress is eudent. . P 
viously stated, tensile transverse stress is a potential source of delamination. As 
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STRESS-Y! ZOOM OF CONTRCT ZONE 




Figure 32. Stress contour: for clamped beam model (In region of contact) 

0(S)-ISO(6)-0(2) case. The resulting smaller contact zone lead to higher contact stresses. 
.-\s seen in the symmetric cases, an increase in the thickness of the layer closest to the 
contact zone resulted in an increase in the magnitude of the ma.ximum tensile transverse 
stress, seen at cross section C. However, the stress at the laminate interface was always 
compressive. 

With regard to bending stresses for these layups, Figures 47 and 48 clearly show the 
nonlinear behavior as the contact layer thickness increases. In addition, the thickness 
of the exterior layer opposite to the contact surface shows similar results as the sym- 
metric cases. As the thickness of this layer decreases, the beam is more susceptible to a 
bending crack that propagates into the interface. 
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Figure 33. Stress distribution: for 0(3)-ISO( 10)-0(3) laminate 
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Figure 34. Stress distribution: t,,, for 0(4)-ISO(8)-0(4) laminate 
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Figure 35. Stress distribution: for 0(5)-ISO(6)-0(5) laminate 
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Figure 36. Stress distribution: for 0(3)-lSO(i0)-0(3) laminate 




Figure 37. Stress distribution: for 0(4)-ISO(8)-0(4) laminate 
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Figure 38. Stress distribution: <r^ for 0(5)-ISO(6)-0(5) laminate 
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Figure 39. Stress distribution: <r, for 0(3)-ISO(10)-0(3) laminate 
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Figure 42. 




Stress distribution: <r, for 0(5)-ISO(6)-0(5) laminate at cross section C 
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Figure 43. Stress distribution: for 0(3)-ISO(6)-0(7) laminate 
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Figure 44. Stress distribution: for 0(8)-1SO(6)-0(2) laminate 
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Figure 45. Stress distribution: <r^ for 0(3)-ISO(6)-0(7) laminate 
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Figure 46. Stress distribution: for 0(8)-ISO(6)-0(2) laminate 
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Figure 47. Stress distribution: o, for 0(3)-ISO(6)-0(7) laminate 



74 




75 



V. CONCLUSIONS AND RECOMMENDATIONS 



This study has developed two methods for approximating contact stresses using the 
augmented Lagrange multiplier method. As illustrated in Part A of Chapter IV, these 
methods accurately approximate the stresses that result from contact between a cylinder 
and plane surface. This study has also illustrated how this approach can be applied to 
understand the beha\ ior of an actual contact problem by examining the response of a 
composite plate to low velocity impact. 

A. RAYLEIGH RITZ APPROACH 

in the process of developing these methods, a number of comments can be made 
regarding the application of the Rayleigh-Ritz method to solving contact stress prob- 



1. The selection of the trial function is an extremely challenging process. If it is de- 
sired to determine the deformation in a contact problem, the proper stress field 
must be first satisfied. Because of this, the selection of possible trial functions is 
limited. For example, when selecting a trial function for the vertical deformation 
of the foundation of Figure 1. a suitable selection is given by the following 
equation: 



This equation exhibits the favorable characteristics of maximum deformation at the 
contact surface and diminishing deformation as the distance from the contact sur- 
face increases. If contact stresses are to be modeled, this trial function is inappro- 
priate. Calculation of is as follows: 



This function exhibits zero strain at the point of contact increasing to maximum 
strain at the lower bounda'ry. 

2. Since the selection of trial functions is difificult, the task is further impeded by 
complicated geometries. Furthermore, selection of a trial function necessitates that 
some knowledge of the deformation field exists. Without a sensible selection of 
trial functions, an accurate approximation is unlikely. 

3. This method assumes the trial function in the form of an infinite series. Solution 
accuracy theoretically should improve with an increased number of terms. How'- 
ever, precautions must be taken to ensure the solution is numerically stable as the 
number of terms increases. Since the strain energy calculations require integration, 
there are choices of trial functions that will increase without bound as the number 
of terms is increased. This problem can be controlled by normalizing dimensions 
or limiting the choice of trial functions. 



lems. 
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4. An increase in accuracy was observed as the number of constraints was increased 
and the distance between consecutive constraints was decreased. It is believed that 
the improved accuracy results from a better definition of the contact surface. 

B. FINITE ELEMENT APPROACH 

The results in Chapter III illustrated that this approach of applying the finite ele- 
ment method to contact stress analysis is effective. A number of comments can be made 
regarding this approach to problem solving. 

1. .As illustrated in the results, this method accurately approximated the isotropic 
roller bearing problem. However, some difficulties were encountered during the 
modeling of the multi-ply composite. In this model, smooth trends of decreasing 
deformations were often interrupted by spurious deformations or groups of defor- 
mations. These interruptions occurred within layers of significantly reduced 
stiffness. It is believed that the optimization routine had difficulty approximating 
the deformations through these layers because of their verx' small contribution to 
strain energy. ,As stated in the results, the contact boundarx conditions were ob- 
tained from the optintization program and applied to a direct finite element pro- 
gram to solve the problem. The above difficulty is recognized as a limitation of this 
approach. 

2. This approach is much more flexible for complicated geometries than the 
Rayleigh- Ritz approach. In addition, detailed knowledge of the deformation field 
is not needed as required by the Rayleigh- Ritz approach. 

3. The application of static condensation is crucial to the successful implementation 
of this method. Every effort should be made to reduce the number of design vari- 
ables to improve optimization efficiency. 

C. COMMENTS ON OPTIMIZATION 

•A number of observations were made regarding the general use of the .Automated 
Design Synthesis System and the specific usage of the augmented Lagrange multiplier 
method. 

1. The global optimum was more likely to be determined when the objective function 
was normalized. 

2. For both the Rayleigh-Ritz approach and the finite element approach, the initial 
choice of design variables had a significant affect on the possibility of obtaining the 
global optimum. For the Rayleigh-Ritz method, initial selections of design vari- 
ables can result in largely dissimilar values of strain energy and external work, the 
two components of the objective function. It was determined that convergence was 
more likely when optimization commenced with these two terms on the same order 
of magnitude. With regard to the finite element approach, sensible choices of the 
initial design variable vector was necessary for convergence to the global optimum. 
This was accomplished by intuitive selection of design variables to model the likely 
deformation. 

3. Solution accuracy can be improved by scaling constraint equations. It has been 
stated that in some circumstances, some constraints change more rapidly than 
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others and can influence the solution excessively while others have little influence 
[Ref. 6: p. 136]. 

4. W^ith regard to the usage of the augmented Lagrange multiplier method, it was 
frequently necessarx’ to 'tune' the optimization algorithm to a specific problem. 
This was done by varying the initial penalty term p and the initial Lagrange multi- 
plier term / . As stated by Vanderplaats, commencing optimization with a small 
value of p should theoretically suffice for most problems [Ref 6: pp. 137-13S]. 
However, it was frequently necessary to select an initial value for p due to conver- 
gence to unrealistic solutions. Similarly , an initial choice of the Lagrange multi- 
plier term can effect the solution. Commencement with a small value is again 
recommended. [Ref 6: p. 141]. This need to tune' the problem is a significant 
drawback to using this optimization method. The ideal way to overcome this lim- 
itation is to first tune the optimization routine using a known solution. With this 
accomplished, this approach can be used for meaningful data collection. 

D. SANDWICH CO.MPOSITE MATERIAL STUDY 

The behavior of sandwich composite materials to low velocity impact loading was 
successfully investigated by the application of the finite element approach. A number 
of observations can be made from examining the results. 

1. The maximum shear stress is relatively insensitive to layer thicknesses. However, 
as the thickness of the contact layer increases, a reduction of the interface shear 
stress is observed. 

2. Tensile transverse normal stresses exist at some cross sections away from the con- 
tact zone. However, this stress is always compressive at the interface. Compressive 
transverse stresses increase in beams with smaller cores due to reduced deflection 
and contact zone size. 

3. As the thickness of the layer closest to the contact zone increases, a nonlinear dis- 
tribution of bending stress within this layer intensifies. This phenomenon is local- 
ized to the region of contact. 

4. .As the thickness of the layer opposite to the contact zone increases, bending crack 
propagation toward the core is less likely due to increased compressive bending 
stresses within the layer. 

E. RECOMMENDATIONS FOR FURTHER STUDY 

The methods developed in this study offer a basis from which additional research 
can grow. ,A reasonable direction is the relaxation of some of the assumptions made in 
Chapter 11 Part B. For e.xample, rela.xation of the rigid roller assumption and the 
frictionless surface assumption w'ould provide challenging research. .Models with com- 
plex geometry could be created. For example, a model of a pin loaded bolt connection 
could be created with rigid or non-rigid pins. Implementation of these changes would 
provide a versatile and highly applicable model for contact stress analysis. 
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